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An analysis technique is presented to provide an essentially explicit 
solution for a system of n simultaneous first-order linear differential equa- 
tions with periodic coefficients. This representation of a periodic variable- 
parameter linear system of arbitrary finite order is chosen for its theoreti- 
cal and. practical advantages over the classical nth. order linear differential 
equation. Emphasis is placed on natural mode solutions of a homogeneous 
set of equations. The characteristic exponents for these solutions are deter- 
mined from a polynomial equation the coefficients of which are linear 
combinations of n — 1 convergent infinite-order determinants. Approximate 
calculation of these determinants is feasible for problems of moderate order. 

I. INTRODUCTION 

Systems of linear ordinary differential equations with periodic coeffi- 
cients are assuming an increasing importance in engineering problems. 
Two applications of present interest are periodically time-variable net- 
works and multimode waveguide with periodic physical distortions. 
Such applications have usually been analyzed by methods appropriate 
to special cases such as the second-order case or by approximate tech- 
niques valid for almost constant-parameter systems. However, pertur- 
bation techniques for almost stationary systems are inadequate for 
careful analysis of large-signal behavior of time-variable networks. 
Similarly, a periodically distorted helix waveguide, for which more than 
two modes must be considered, 1 should be described by a differential 
system of order greater than two. These examples illustrate the impor- 
tance of a technique for obtaining essentially explicit solutions of 
periodic variable-parameter linear systems. Solutions in terms of char- 
acteristic exponents are known to exist for systems of linear differential 
equations with periodic coefficients. 2 However, the methods usually 
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employed for solving such systems, such as power-series techniques, 
iterative processes, and incremental numerical solution methods, fail 
to provide a system response description valid for all values of the 
independent variable (time, distance, etc.). 

The analysis method to be presented below provides an essentially 
explicit solution for periodic variable-parameter linear systems of arbi- 
trary finite order. The solution describes the system behavior for all 
values of the independent variable. Emphasis will be placed on obtain- 
ing a set of basis functions for a homogeneous system, since the solution 
in the inhomogeneous case can be obtained from the basis functions. 
As shown by Darlington, 3 these functions may be regarded as analogues 
of partial fractions in fixed network theory. 

II. FORMULATION OF DIFFERENTIAL SYSTEM 

In this discussion the system of equations to be solved will be repre- 
sented by the vector differential equation 

F'(t) = B(t)F{t) (1) 

where F{t) and B(t) are nth-order column and square matrices, re- 
spectively, and the prime denotes differentiation with respect to the 
independent variable t. It is supposed that the elements of B(t) are 
known functions of t with a common period of unity, i.e., 

B(t) = B(t+ 1). (2) 

The formulation of this problem in (1) is chosen not only for its 
elegance, but also because of its practical advantages. As indicated by 
Kinariwala 4 these include the ability to write such an equation directly 
from a time- variable network, the fact that the eigenvalues of B{t) are 
natural frequencies for stationary networks, and the convenience of (1) 
in obtaining the quadratic forms representing stored energy and dissi- 
pated power in stationary or nonstationary cases. These advantages have 
their translated versions in other physical problems, including multi- 
mode waveguide problems. Moreover, an equation such as (1) is 
easily obtained from an nth-order linear differential equation, but the 
transformation from (1) to such an equation can be quite difficult (or 
analytically inconvenient). 4 Thus, (1) represents a well-founded be- 
ginning for the analysis of variable-parameter problems of practical or 
theoretical interest. 

III. FORM OF SOLUTIONS 

The form of solutions of (1) is well known; 2 pertinent properties of 
such solutions will be reviewed here briefly. If B(t) is piecewise con- 
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tinuous (1) has the unique solution 

F(t) = X(t)F(0) (3) 

where X(t) is the unique nonsingular square matrix satisfying 

X' = BX 

(4) 

X(0) = I = diag{l}. 

When B(t.) satisfies (2), X(l) may be written as 

X{t) = J(t) e Kt (5) 

where 

J(t) = J(t+l) (6) 

and 

e K = X(l). (7) 

For convenience it will be assumed here that the eigenvalues of K 
are distinct, or at least that K can be diagonalized; thus, a constant 
nonsingular matrix P exists so that 

K = PMP' 1 (8) 

where 

M = diag {n,\ (9) 

and the constants m are the eigenvalues of K. The matrix exponential 
function in (5) may be similarly diagonalized, so that the solution (3) 
may be constructed in the form 

F(t) = J(t)P[dmg { e ""}]P-^(0). (10) 

By establishing the special initial conditions 



F,(0) = P 



i th row 



(11) 



the corresponding unique solution 

Fi(t) = e't'JWFm, (12) 

is obtained from (10). Thus, by proper choice of initial conditions a 
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set of n solutions of the form 

F(t) = e" l Q(t) (13) 

where n is a scalar constant and Q(t) is a column matrix with period 
unity, have been shown to exist. 

The n solutions in the form (12) or (13) represent natural modes of 
the periodic system described by (1) and (2). If the n values of m are 
distinct the corresponding n solutions are certainly independent and 
form a set of basis solutions of ( 1 ) . Any other solution of ( 1 ) comprises 
a linear combination of solutions like (12) or (13). Moreover, as Dar- 
lington 3 has pointed out, these natural-mode solutions are essentially 
unique because of their simple form. Hence the natural modes given by 
(12) represent a complete and essentially unique description of the 
natural behavior of the periodic system. The eigenvalues m , frequently 
referred to as characteristic exponents, play a role analogous to response 
poles or natural frequencies of stationary systems. The strength of each 
natural mode in the homogeneous case is determined by the initial 
conditions and the constant matrix P. Moreover, the natural-mode 
solutions allow a complete solution to be calculated in the inhomogene- 
ous case. Thus, the determination of n corresponding solutions for /* 
and Q(t) in (13) is central to the problems associated with (1) and (2). 

The object of the present treatment is to indicate a technique for 
determining the characteristic exponents n, as well as the corresponding 
matrices Q(t) if desired. Primary attention is given in finding the char- 
acteristic exponents n because of their practical importance and because 
the solution for Q(t) is not greatly difficult in principle if the appropriate 
characteristic exponent is known. Solutions for Q(t) are mentioned in 
Appendices A and B. 

The method to be discussed resembles the technique used by Hill 5,6 in 
solving the second-order equation 

x"(t) + A(t)x(t) = (14) 

where A(t) is periodic. It will be shown that the characteristic exponents 
may be determined from roots of either a transcendental or polynomial 
equation in which certain infinite-order determinants enter as parame- 
ters. A technique similar to Hill's was employed by H. von Koch in the 
last century to provide an explicit solution in terms of infinite-order 
determinants for a general nth-order linear differential equation with 
periodic coefficients. This technique is carefully discussed by Forsyth 
and Riesz, 8 who also give references to von Koch's original papers. Thus, 
the method presented here, although developed independently, does 
not solve an unsolved mathematical problem when applied to a periodic 
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variable-parameter system described by an nth-order linear differential 
equation. It does, however, solve the stated problem in a way that 
appears to have several advantages, mostly associated with its formula- 
tion as a system of n simultaneous first-order linear equations. These 
advantages, already mentioned in Section II, seem likely to make the 
present solution technique more useful in the analysis and synthesis 
of periodic variable-parameter systems than one based entirely on the 
classical nth-order linear differential equation. 

IV. INTEGRAL FORM FOR THE PERIODIC SYSTEM 

The analysis of the periodic system begins by multiplying both mem- 
bers of (1) by e~ at , where a is an arbitrary constant, and adding and 
subtracting aFc~"' to yield, whenever F' exists, 

(/<>-")' + aFe~ at = BFe~ at . (15) 

Integration of (15) results in the integral equation 

Fc-"' + a f Fe~ at fit = J BFe~ at dt + C (16) 

where C is a constant. Any solution of (15) is also a solution of (Hi); 
thus, let /'' be a solution given by ( 13) and let 

a = /x +./2xA- (17) 

where /.• is an arbitrary integer. Equation (16) becomes 

Qc-^' + ( M + ./'27T/.0 f Qe'^'dt = f BQe~ j ' iTkl dt + C. (18) 

If (18) is evaluated at / = and t = 1, and the results subtracted, the 
first term in (18) makes no contribution, being periodic. Hence, (18) 
implies 



( M + j2 v k) f Qe-* wkt dt= f BQc- j2 ' kl 



dt (19) 



for all integers /,-. It will be seen below that this integral equation suffices 
to determine n and Q(t), which are essentially eigenvalues and eigen- 
f unctions. 

V. MATRIX DIFFERENCE EQUATION 

To make use of (19) in finding solutions of (1) it will be assumed that 
the given matrix B(t) and the solution matrix Q(t) may be expanded in 
the Fourier series 
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B{t) = E B p e iUpl (20) 



p= — oo 

and 



Q(t) = E Q,«*" 



(21) 



where matrices i? p are square matrices and Q p are column matrices. 
Requirements on the asymptotic behavior of the elements of matrices 
B p and Q p for large values of | p | will be discussed in Appendix A in 
relation to convergence of certain infinite-order determinants. The 
Fourier series for the matrix product BQ may be written as 

BQ= E (BQ) p e i2 * pt (22) 

p=— 00 

in which the column matrices (BQ) P are given by the convolution 

X 

(BQ) p = E fip-rQr- (23) 

r=— oo 

Except for a factor of 2x the integrals in (19) express the Fourier 
coefficients of Q and BQ. Thus, if Q and BQ possess Fourier scries (19) 
is equivalent to the infinite set of linear equations 

( M + J2**)Q» = (BQ), (24) 

or 

CO 

U+ja*)0i = E A-rQr, (25) 

when? /; assumes all integral values. Equation (25) might be regarded 
as a matrix difference equation for Q p ; however it is more convenient 
here to consider (25) as defining an eigenvalue problem for an infinite 
matrix. In terms of Kronecker's 8, (25) is 

00 

= E Wk-r - Mm + i27r/c)/]Q r (26) 

T=— 00 

where / is the nth-order unit matrix. The expanded form of (26) is 
shown in the following infinite-order matrix equation, in which the 
first matrix is partitioned into n X n size blocks and the second into 
n X 1 size blocks. The "origins" of the matrices fall at (B — nl) and 
Qo. 
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Bo - (m - ;4x)7 B-i 5_ 2 

£, Bo - 0* - J'27r)/ B_, B-5 

£ 2 £, Bo -ill B-i B-» 

B s Bi Bo - (v + j2tt)I B-i 

B-i B x Bo-Ui+j**)! 



Q-i 







Q-i 







Qo 


= 





Qi 







Q, 








For convenience it will be assumed that B is a triangular matrix so 
that its eigenvalues appear explicitly as main diagonal elements. To 
show that a constant linear transformation of the dependent variable 
F can always produce this property, let 



B(t) = B + A(t) 
where A(t) has a zero mean, and let 

X = PF 



(28) 



(29) 



where P is a nonsingular matrix of constants. Then ( 1 ) is transformed to 

X' = {PB,P~ l + PAP-')X. (30) 

This equation has the same form as (1), but the constant term in its 
coefficient matrix is the matrix PBJ*~~ l derived from B by a similarity 
transformation. It is well known that a square matrix is reducible by a 
similarity transformation to the classical canonical form having eigen- 
values on the main diagonal and possibly nonzero constants in some 
positions of the next higher diagonal. 9 (These constants cannot appear 
if Bo has distinct eigenvalues; hence, B can often be assumed to be 
diagonal.) The matrix B can also be reduced to a triangular form by a 
similarity transformation in which P is a unitary matrix. 1 . This re- 
duction, which is always possible, may sometimes have advantages in 
studying energy functions or related quadratic forms. Thus, by either 
technique B u can be reduced to triangular form. It will be assumed that 
such a transformation has been effected in obtaining ( 1 ) . 

VI. CONVERGENT LINEAR EQUATIONS AND INFINITE-ORDER DETERMINANT 

To produce convergence of the determinant of coefficients of the 
infinite set of homogeneous equations denning Q, Equations (26) or 
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(27), divide each elementary row of the coefficient matrix in Equation 
(27) by its main-diagonal elements. When B is diagonal this process is 
identical with multiplication of each equation in Equation (26) or 
each matrix row in the square matrix of Equation (27) by the diagonal 
matrix [B — (n + j2vk)iy l . In general, let the matrix A be defined by 



A = diag {X p { 



(31 



where X p represent the n eigenvalues of B . Then the set of equations 
with convergent determinant may be written as 

[Mtr)[Qr] = (32) 

where submatrices Mkr are given by 

Mkk = [A - (n + j2Tk)I]- 1 [B - (> + j2wk)I] (33) 

and 

M kr = [A - (m + j2irk)IY l B k . r , k * r. (34) 

When Bo is diagonal, Mkk reduces to the unit matrix. Thus, a typical 
determinant d[M kr ] of (32) may be illustrated for n = 2 by the following 
scheme : 



d(») = 




ai 


1 

6i 


Xi — n 
Cl 


Xi — M 


X 2 — m 
a* 


x 2 - M 
b 2 


Xi—j2ir — n 
Ci 


\i—j2ir — n 



et-i 


6_i a_2 6_2 


M+j2ir- 


"M Xi+j'2jt-m Xi+J2tt-m Xi+J^ir— m 


C_l 


d_l C_2 f/_2 


\2+j2w- 


-/i X 2 +j2jt— m X2+j2jt— ju Xj+j2x— m 




a_! 6_i 




Xi — /z Xi — p 





C_i d_i 



Oi 



&, 



X 2 — M 

1 



X 2 - M 




X'—firr — ft \-i—j2ir — fi X 2 — j2nr— fi \»—j2w — fi 

The notation used in (35) is 

p ^ 



o 







a p 


*S 




#, 




A 


d P _ 


> 






X, 







Bo = 










Lo 


xj 



(35) 

(36) 

(37) 



SYSTEMS OF LINEAR DIFFERENTIAL EQUATIONS 1283 

The determinant of the infinite-order matrix [M kr ] of (32), illustrated 
by (35) for n = 2, will be denoted by d(p) to show its functional de- 
pendence on the argument ft. The function d(n) is actually a determi- 
nant of infinite order. If this determinant converges it represents a 
function of n which must vanish in order to obtain nontrivial solutions 
for Q r in (32). Requirements necessary for the convergence of d(ji) are 
discussed in Appendix A, where it will be shown that d(n) converges 
for a large class of problems. Hence the basic equation 

d(/0 = (38) 

defines the characteristic exponents of the differential system (1). 

VII. FUNCTIONAL EXPRESSIONS FOR THE CHARACTERISTIC DETERMINANT 

Equation (38) taken alone is rather unwieldy, involving as it does 
the equation to zero of an infinite-order determinant whose elements 
are functions of n. However, it will now be shown that expressions for 
d(n) in terms of elementary functions may be written to allow a simple 
solution of (38). 

The determinant d(n) is shown in Appendix A to converge for all 
values of n except those for which the denominators of rows of d(n) 
vanish. Multiplication of one row of an infinite-order determinant by any 
scalar is equivalent to multiplication of the determinant by the same 
scalar. Similarly multiplication of any row of d(fi) by its corresponding 
denominator X ; , — (m + j2vk) produces a determinant convergent at 
\ p = p + f2wk, so that each row of d(n) introduces exactly one pole 
in d(n). Moreover, d(n) is periodic in n with period $2ic, since replacing 
H by n + j2ir only shifts the origin of the infinite-order determinant. 
Evidently d(n) has simple poles at 

M = a p + j2irq, p - 1,2, • • •, n q integral. (39) 

It will be assumed for the moment that these poles are distinct; this 
restriction may be relaxed slightly, as shown in Appendix B. Finally, 
as n approaches infinity along any radial line in the complex n plane 
except a vertical line, the off-diagonal elements in d(n) tend toward 
zero, or briefly 

rf(x) = 1. (40) 

The periodicity of d(n) implies that the residue of d(n) at any of the 
poles in (39) is independent of the particular integer q. Thus, a formal 
expansion of d(n) in partial fractions is 

d(n) = K« + ± E -*g .. . (41) 

p =i q — * n — \ p — jZirq 
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According to the Mittag-Leffler theorem 11 this expansion defines the 
function 



<*<M) = K« + \ £ K P coth (^- P ) • 



(42) 



A relation fixing K x may be derived from (40) by noting that as n 
approaches infinity along any nonvertical radial line 

Hm coth (j-^) - 1 (43) 

so that 

L-l-iEZ,. (44) 

To compute the residues if p the well-known rule 

K p = lim (it - \ p )d(ix) = [(/i - X p )d(/t)L-x p (45) 

M-»X p 

is employed. The procedure is simply to multiply every element in the 
row of d(n) containing X„ — n (in the denominators) by the factor 
(n — X p ) and to evaluate the resulting determinant. For example, in 
the case of n — 2 used above for illustration, the row of d(\i) containing 
Xi — \i in the denominators is replaced by 

• • • — a 2 — 62 — ai — hi — a_i — 6_i — a_2 — &-2 ■ ■ ■ (46) 

and the resulting determinant evaluated at m = Xi . Reasonably accurate 
and efficient methods for computing K p from such a determinant can 
be programmed just as for Hill's determinant in the second-order case. 
Such a technique is discussed briefly in Appendix C. 

It is well known that the solution of Hill's equation generally requires 
the evaluation of only one infinite-order determinant, while the solution 
of a second-order problem using (38) and (42) appears to require the 
evaluation of two determinants. Actually it will be shown that only 
n — 1 determinants need be calculated for an nth-order system of equa- 
tions. In addition (42) may be simplified because of the relation among 
the residues K p to be demonstrated below. 

To examine the poles and zeros of d(n) it is convenient to consider 
the complex fi plane divided into horizontal strips of width 2w. The 
poles of d(n) fall at X p + j2rq. Although the eigenvalues X p may lie 
in any of these strips, values of q always exist to give one pole in the 
fundamental strip ^ Im m < 2x representative of each X p . Hence 
d(n) has exactly n poles in each strip. It will be seen shortly that d(n) 
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also has n zeros in each strip so that a pole-zero constellation for d(n) 
might be illustrated by Fig. 1. 

The desired relation among the residues K p is obtained by noting that 



/ 



d{n) dn = £ K t 

abed p=l 



(47) 



where the integral is taken around the rectangular contour abed shown 
in Fig. 1 (or a congruent rectangle vertically displaced if a pole happens 
to fall at Im p. = 0). The periodicity of d{n) insures that the contri- 
butions to the integral from the horizontal sides ab and cd will cancel. 
The vertical sides be and da are supposed to be displaced from the 
origin far enough to include all n poles in the rectangle so that (47) is 
valid. As their displacement approaches infinity the value of d(n) ap- 
proaches unity and the contributions from the vertical sides tend to 
cancel. Thus (47) implies 

(48) 



E K p = 0. 
This relation shows that (44) and (42) may be simplified to 

K« = 1 

and 



(49) 



din) = 1 + I E K p coth (j—^) ■ (50) 

It also allows one residue to be computed from the other n — 1, al- 

Im/j. 
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Fig. 1 — Pole-zero constellation. 
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though all n residues might be computed in practice and (48) used as a 
check for numerical accuracy. 

Equation (50) expresses the characteristic determinant d(/x) in 
terms of the eigenvalues X p of the stationary part of the system and the 
residue determinants K p . The characteristic exponents n are thus by 
(38) and (50) the roots of the trigonometric equation 

= 1 + |g^coth(^^). (51) 

This trigonometric equation represents an explicit solution of the prob- 
lem of rinding characteristic exponents for an nth-order periodic system. 

It is evident from (50) as well as from Fig. I and the periodicity of the 
function e* 1 that the substitution 

z = e* (52) 

reduces (50) to a rational function in z. Zeros and infinities of z do not 
introduce superfluous poles or other singularities in this function be- 
cause of (43). Thus, the poles and zeros of this rational function are 
mapped by (52) into the poles and zeros of d(n) shown in Fig. 1. Any 
strip of vertical width 2ir in the \i plane is mapped by (52) into the entire 
z plane so that the rational function of z has n poles in the z plane. The 
number of zeros of the rational function is also necessarily n. Hence, 
d(fi) has precisely n zeros in any horizontal strip of width 2x in the n 
plane. 

Because of the existence of well-developed techniques for polynomial 
manipulation, such as approximate solution methods, interpolation 
formulae, and stability criteria, it is practically convenient to utilize 
(50) and (51) in rational form. Accordingly let z be defined by (52) and 
x P by 

x p = A, V = 1,2, •••,n, (53) 

so that d(n) is transformed to 

D(z) -l+^Jf, (i±^A = tfdog 2 ). (54 ) 

2 p =i \z - x p / 

Further, let 

g{z) = fl (z - x P ) (55) 

be a characteristic polynomial defining the eigenvalues of the stationary 
part of B(t). (This "characteristic polynomial" differs from the con- 
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volitional one in that its roots are c Xp rather than X p .) Equation (51), 
the characteristic equation, then becomes 

= /(z) + g(z) (50) 

where 



/(«) =lZ K p {z+.r P )fl(z -•«•„). 



y=l 
q^p 



(57) 



These equations demonstrate that the characteristic polynomial for 
the periodic system is obtained by adding a certain interpolating 
polynomial to the characteristic polynomial of the stationary part of 
the system. The behavior of the interpolating polynomial is prescribed at 
the roots of the stationary part. 
The interpolating polynomial /(z) has the n assigned values 

f(x p ) = KpXp II (* P ~ * a ); (58) 

'7=1 

•I^P 

because of the relation (48) among the residues K p the polynomial 
(57) is identical with the Lagrangian interpolating polynomial 

/(-') = E K„.v P tl (z - x g ). (59) 

p=i </=i 

Evidently the interpolating polynomial /( z) is the unique polynomial 
of degree n - 1 having the assigned values (58). Thus/(z) + g(z), the 
characteristic polynomial of the periodic system, is the unique monic 
polynomial of degree n having the n assigned values given in (58). 
This point of view may give some insight into stability questions. 
For example, the classical criteria of Routh and Hurwitz, and other 
results on bounds of zeros of sums of polynomials may be useful here. 

If all the residues K p vanish, as in the stationary case, the limiting 
values of the characteristic exponents obtained from (51) and (50) 
are u = X p , p = 1 ,2, • • • , n. In cases of small variations where all 
| K p | are small the characteristic exponents differ very little from the 
eigenvalues of the stationary part. Asymptotically they may be calcu- 
lated from any of the approximate equations 

0%l + iA>oth^-^) (60) 

z«x p (] - K p ) (01) 
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or 



X — K v . 



(62) 



Although perturbation type solutions such as (62) probably are more 
easily calculated by less complicated techniques, characteristic expo- 
nents obtained from (61) or (62) may be useful as starting values for 
solving (51) or (56) by numerical methods. 

VIII. EXAMPLE 

The following example illustrates the technique for finding charac- 
teristic exponents. A second-order case is chosen for convenience 
because some digital computer programs needed for the efficient evalua- 
tion of the residue determinants are not yet available. However, higher- 
order examples are not different in principle nor will they require in- 
ordinately longer computations. 

The Mathieu equation 



dy 



+ (3 - 4 cos 2z)y = 



has the solution 



with 



= e'* J2 C2r+l e 



i(2r+l)z 



= ±0.57943224- • 
In vector form this equation is equivalent to 





r = 7T 



i 



-3 + 4 cos 2vt 



with the identifications 
Y = 



£]■ 



2/i = y, 



z = irt. 



Diagonalization by the transformation PY = F where 

i r jV3 i" 



p = - 



#v3L-iV3 i 



(03) 

(64) 
(65) 
(66) 

(07) 

(08) 



yields (1), with 
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m -m 



•-* , +in/5 



1 








_0 -1_ 






+ JV3 


|~1 
1 


-ll 
-1_ 



(69) 



»j2»* 



The determinant d(n) has the form shown by (35), and the residue at 
/x = ,/""v/3 is approximately K = j'0. 562096, a result obtained from a 
42nd-order approximant. From (56) and (57) the characteristic ex- 
ponents are solutions of 

cosh n = cos 7r\/3 — jK sin 7r\/3, (70) 

which yields, for K « J0.562096 

M « ±j'0.42059ir. (71) 

Corresponding correct values of p. from (64) and (65) are ±j0.420577r. 
A somewhat longer computation would be required to produce a re- 
sult accurate enough for certain purposes. Such a computation was not 
employed here because a more fundamentally sound computing tech- 
nique for band-limited periodic variations as in (69) would exploit 
the form of the residue determinant and its large number of zero ele- 
ments. Specifically, it is possible to program a determinant evaluation 
technique for such cases so that the computation time is asymptotically 
proportional to the order of the truncated determinant rather than to 
its cube. This possibility is discussed further in Appendix C. 

IX. CONCLUSIONS 

A method has been developed for analysis and calculation of solutions 
of nth-order linear periodic differential systems. The system description 
employed is a set of n simultaneous first-order linear differential equa- 
tions. The method allows the determination of characteristic exponents 
from polynomial equations the coefficients of which are linear combina- 
tions of n — 1 convergent infinite-order determinants. Approximate 
computation of the determinants is feasible for problems of finite order. 
In addition to characteristic exponents the complete solutions may also 
be computed if desired. 

APPENDIX A 

Convergence 

The validity of the analysis presented here depends upon the con- 
vergence of the infinite processes employed. It must be shown that the 
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determinant d(n) and the Fourier series for Q(t) are convergent if the 
coefficient matrix B(t) is suitably restricted. For this purpose (32) may 
be written as the infinite set of scalar equations 

■x. 

Xi + Z otto = (72) 

where a,:; and xj are scalars, and the equations hold for all integral i. 
The coefficients a,, actually are elements of the submatrices M kr , and x t 
elements of submatrices Q r in (32). The determinant of coefficients of the 
scalar equations is 

d(/x) = | hi + an | . (73) 

According to a theorem of St. Bobr 13 this determinant is absolutely 
convergent if 

E \au\ (74) 



and 

00 

S 



SI c |P/(P-D 
I U 'iJ I 



_ I?*! 



(75) 



converge for some value of p in the interval 1 < p ^ 2. (For p = 2, 
the case used here, the theorem was given by von Koch.) The expression 
in (74) obviously converges to zero since all a,-,: in (72) are zero. Let the 
elements of the given matrix B(t) be square integrable functions. Then 
Parseval's relation applies and the Fourier series coefficients for the 
matrix elements are surely square summable. Hence, the inside sum 
in (75) converges for p = 2. The outside sum also converges for 
p = 2, since its general term is asymptotically proportional to i~ 
for large \i\ (as (33) and (34) indicate by their dependence on A;). 
Of course, an exception occurs for values of n given by Equation (39). 
The determinant d(/i) is singular at these points, but the convergence 
of the residue determinants K p for simple poles is assured by St. Bobr's 
theorem. Thus, d(n) converges absolutely and uniformly except for n 
arbitrarily near X p + j2irq and has poles at these values of \x. 

Since the determinant d{n) has zeros at any of the n characteristic 
values of n within the strip ^ Im,\i < 2ir, the deletion of the zeroth 
equation (i = 0) from (72) and the transposition of an&o in each equa- 
tion produces a nonzero determinant of coefficients in the equations 

1 Xi + E dijXj = — aupct = iji, i 7* 0. (76) 
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These equations with a, 7 evaluated at a characteristic value of n have 
as solutions the scalar quantities .r, needed to produce the matrices 
Q p and thus the matrix Q(t). That meaningful solutions to (76) exist 
for an arbitrary constant .To is shown by a theorem of L. W. Cohen. 
This theorem (paraphrased) states that if (75) converges for the 
coefficients in (76), if the (convergent) determinant of (76) does not 
vanish, and if 

00 

i=—x 

converges, then the solutions exist, may be obtained by Cramer's 
rule (with infinite-order minor determinants), and have the property 
that 



£ I x i 



p 



converges. Thus, if the elements of the given matrix B(t) are square 
integrable functions, the coefficients a l0 are surely square summable, 
and the resulting trigonometric series for elements of Q(t) have square 
summable coefficients. The Riesz-Fischer theorem 1 then states that the 
elements of Q(0 are square integrable functions with Fourier coefficients 
given by the elements of Q p and that the Fourier series for Q(t) con- 
verges to Q{t) in the mean. (Consequently there exists a sequence of 
partial sums of the Fourier series converging to Q(t) "almost every- 
where.") 

In a more restricted case which might have more practical importance, 
it may be shown that if B(t) is continuous so that the elements of B p 
are 0(]/p 2 ), the solution matrix Q(t) has the same property. Of course, 
the Fourier series for Q(t) converges absolutely and uniformly in this 
case. This convergence condition and the more general one above demon- 
strate that the analysis technique is valid for a wide class of problems. 

APPENDIX 13 

Multiple Poles of the Characteristic Determinant 

If the matrix B , the stationary part of the coefficient matrix B(t), 
has repeated eigenvalues, or if any of its eigenvalues differ by integral 
multiples of j2ir, some denominators of rows of d(n) are identical. In 
this case d(n) has multiple poles, and the necessary analytical and 
computational procedures become more complicated. It is possible to 
treat the case of a single second-order pole of d(jx) by evaluation of n — 1 
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determinants as before, but greater multiplicities require considerably 
more extensive calculations. 

When d(n) has an ra-fold pole at /j = Xi the partial fractions expan- 
sion of d(n) must contain the corresponding principal part of d(n). 
The coefficients in the principal part involve derivatives of (n — Xi) m 
d(n) evaluated at /x = Xi . These derivatives are more difficult to 
compute than the residue determinants of the simple case because they 
are linear combinations of most of the first minors of d(n). The com- 
putation of such minors (not necessarily by direct methods) is also 
required if Q(t) is to be determined (even when d(n) has only simple 
poles). Appendix A shows this computation to be theoretically possible; 
it is equivalent to the inversion of a set of equations like (76). Never- 
theless, the computation effort would be considerably greater than that 
required for computation of characteristic exponents when d(n) has 
only simple poles. 

When d(n) has a single second-order pole, (48) may be utilized to 
make possible the calculation of characteristic exponents. It is conven- 
ient here to use the rational form of (54) for the infinite-order determi- 
nant d(n). Let the repeated roots of B be identified with Xi , X 2 and 
X\ . x 2 respectively. Define ai and ai by 



a x = Ki(xi — x 2 ) 
a 2 = K 2 (x 2 — x^ 
and allow X\ to approach x 2 . Substitution of (77) in (54) yields 



(77) 



2 \xi — xj (z — Xi)(z — x 2 ) 

2 P =3 V - x p/ 



(78) 



+ 



2(Z - Ei) (2 — Xi) 

Equation (48) may be written as 



ai — a 2 



Xl — X 2 p=3 

so Hint 



+ £ K p = (79) 



where 



L + Z lim K p = (80) 

p=3 H-+ZO 



L = lim fa °*) = lim (Kt + K 2 ) (81) 

x l~* x i \ X l •&*/ x \~* x 2 
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is a finite limit. Evidently, as X\ approaches x 2 , 
L (z + x 2 



a? 



+ \£ K '(£ir)- (82) 



HmZ)(s) = l+ 5 (^j + (z _ j2)! . 2; _ 

The zeros of this limiting form of D(z) correspond to the characteristic 
exponents in this case. The parameter a 2 may be determined by factor- 
ing 1/(X2 — Xi) from the appropriate row of K 2 and computing the 
resulting determinant A, since 



«2 



= A-lim 



lim (r^f ) 

l-x 2 \Al A2/ 



= A. 



(83) 



The parameter L required in (82) may be computed from (80). Cases 
where two poles of d(n) are almost coincident may be treated in a 
similar fashion, except that no limits are involved. 



appendix c 



Approximate Computation of Residue Determinants 

In practical cases where the number of terms in the Fourier series for 
B(t) is limited, truncated approximants to the residue determinants may 
be evaluated by techniques that exploit the special form of these de- 
terminants. The form of a truncated residue determinant is illustrated 
by the scheme in Fig. 2, in which all elements outside of the shaded 
region arc zero. Except for one submatrix near the center of the array 
the principal diagonal blocks represent nonsingular triangular sub- 




Fig. 2 — Form of truncated residue determinant. 
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matrices. To demonstrate the feasibility of computing truncated residue 
determinants of large order it will be shown that the computation time 
required for a reduction to triangular form is much smaller than for a 
general determinant of the same order. (The computation time required 
for a general determinant is asymptotically proportional to the cube of 
its order.) 

To evaluate the determinant in Fig. 2 let zeros be produced below 
submatrix 1 by elementary operations with the rows passing through 1 . 
Similar operations to produce zeros below 2 do not disturb the zeros 
already produced. Such operations may be continued in the usual manner 
to produce zeros below 3, 4, etc., until a triangular array of submatrices 
is realized. The number of arithmetic operations necessary in each step 
of zero production is essentially dependent only upon the order of the 
original system of equations and the number of terms in B{t). Observa- 
tion of Fig. 2 shows that the number of zero-producing steps for a 
truncated determinant of large order is asymptotically proportional to 
the order of the determinant. Thus, the computation time required for a 
reduction to triangular form is also asymptotically proportional to the 
order of the truncated determinant to be evaluated. 
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